% This file can generate Figure 3, Figure 4, and Table 3 in the paper. 
% Also, the file is the one generating the results needed to reproduce
% Figure 7. In order to obtain the latter, this file needs to be run
% several times, one for each measure of output (option 'use_data' below);
% each time the file is being executed, it will automatically save the
% results in a .xlsx file that later will be read in the file named
% 'FIG_7_Robustness_Data_Output.m'. 

close all; clear; clc;

global P2 lambda 

%% Options for the specific results to generate.

use_data    =   1;              % = 1 if baseline data
                                % = 2 if output is Industrial Production Index 
                                % = 3 if output is hours per capita nonfarm
                                % = 4 if output is unemployment rate.

csv_file  	=   1;              % = 1 if you want to import from CSV.
                                % = 0 if you want to import from XLSX.
                                
graph       =   1;              % = 1 to generate a plot. 

table       =   1;              % = 1 to generate the first two columns of 
                                %     Table 3.
                                % = 2 to generate the pointwise test
                                %     (third column in Table 3).
                                % = 3 to generate the pathwise test (last 
                                %     column in Table 3)

%% Importing and Preparing the data.

if csv_file == 1
    data            =   csvread('Data_File_Baseline.csv');
    % UM Expectations Median:     
    pilev_data      =   csvread('Data_File_Inflation.csv');
    pilev           =   pilev_data(:,6);     
else
    data            =   xlsread('Data_File.xlsx','Baseline');
    % UM Expectations Median:
    pilev           =   xlsread('Data_File.xlsx','Inflation','F2:F281');
end
dz          =   data(:,3);
dy          =   data(:,4);
zlb         =   data(:,5);
rec         =   data(:,6);
zlev        =   data(:,7);
time        =   data(:,9);
[T,n]       =   size(dy);

if use_data == 1             % GDP
    ylev        =   data(:,8);      
elseif use_data ==  2        % IPI      
    ylev_data   =   csvread('Data_File_Output.csv');
    ylev        =   log(ylev_data(:,2));  
elseif use_data ==  3         % Hours
    ylev_data   =   csvread('Data_File_Output.csv');
    ylev        =   (ylev_data(:,6));  
elseif use_data ==  4          % Unemployment
    ylev_data   =   csvread('Data_File_Output.csv');
    ylev        =   log(ylev_data(:,8));  
end

% Defining the sample:
st          =   148;                        % = 1 for full sample 
                                            % = 148 for 1984Q1 
                                            % = 52 for 1960Q1. 
%% Options for local projections.

P           =   2;      % Number of lags for endogenous variables.
H           =   9;      % Maximum number of leads.
h1          =   1;      % For splines.        
lambda      =   100;    % Value of penalty of SLP.
r           =   3;      % r-1 = order of the limit polynomial for the SLP.

%% Creating variables.

% LHS:
ylev2       =   ylev(st:end);
% Productivity:
zlev2       =   zlev(st:end);
% Inflation:
pilev2      =   pilev(st:end);
% Dummy for ZLB:
zlb2        =   zlb(st:end);

% Creating W matrix: 
ff          =   [ylev2 zlev2 pilev2];
w           =   lagmatrix(ff,1:P) ;     % Remove period t observations.

% Creating lagged variables interacted with period-t ZLB dummy:
[TY,YT]     =   size(w);
zlb3        =   zeros(TY,YT);
for jx = 1:YT
    zlb3(:,jx)  =   zlb2;
end
www         =   zlb3.*w; 
zlb4        =   1-zlb3;
wwww    	=   zlb4.*w;

%% Creating the IRFs.

% LHS variable:    
y                       =   ylev2(1+P:end,:); 

% Local projection.
        % Standard case.               
w                       =   [zlb2(1+P:end) zlb2(1+P:end).*zlev2(1+P:end),...
                            wwww(1+P:end,:), www(1+P:end,:)];
x                       =   (1-zlb2(1+P:end)).*zlev2(1+P:end,:);  % Shock variable.    
lp                      =   locproj_var(y,x,w,h1,H,'reg'); 
if table == 3
    lp                      =   locproj_var_sur(y,x,w,h1,H,'reg'); 
end
% Extract the IRF from local projection: 
irf                     =   lp.IR(2:H+1);
Xp                      =   size(lp.X,2)/H;
% Construct standard errors:
u                       =   lp.Y - lp.X*lp.theta;
[nwse_lp,vc_lp]         =   NeweyWest(u,lp.X,0,0);
    % ZLB case.
w                       =   [zlb2(1+P:end) (1-zlb2(1+P:end)).*zlev2(1+P:end),...
                            wwww(1+P:end,:), www(1+P:end,:)];
x                       =   zlb2(1+P:end).*zlev2(1+P:end,:);   
lp_zlb                  =   locproj_var(y,x,w,h1,H,'reg'); 
if table == 3
    lp_zlb                  =   locproj_var_sur(y,x,w,h1,H,'reg'); 
end
% Extract the IRF from local projection: 
irf_zlb                 =   lp_zlb.IR(2:H+1); 
% Construct standard errors:
u                       =   lp_zlb.Y - lp_zlb.X*lp_zlb.theta;
[nwse_lp_zlb,vc_lp_zlb] =   NeweyWest(u,lp_zlb.X,0,0);

% Smooth local projection.
        % Standard case.  
w                       =   [zlb2(1+P:end) zlb2(1+P:end).*zlev2(1+P:end),...
                            wwww(1+P:end,:), www(1+P:end,:)];
x                       =   (1-zlb2(1+P:end)).*zlev2(1+P:end,:); 
slp                     =   locproj_var(y,x,w,h1,H,'smooth',r,lambda);
if table == 3
    slp                     =   locproj_var_sur(y,x,w,h1,H,'smooth',r,lambda);
end
% Extract the IRF from SLP: 
irf_slp                 =   slp.IR(2:H+1);
P2                      =   slp.P2;
% Construct standard errors:
u                       =   slp.Y - slp.X*slp.theta;
[nwse_slp,vc_slp,Q]     =   NeweyWest_slp(u,slp.X,0,0);
V2                      =   slp.B*vc_slp(1:slp.K,1:slp.K)*slp.B';
se_slp                  =   sqrt(diag(V2));
        % ZLB case.
w                       =   [zlb2(1+P:end) (1-zlb2(1+P:end)).*zlev2(1+P:end),...
                            wwww(1+P:end,:), www(1+P:end,:)];
x                       =   zlb2(1+P:end).*zlev2(1+P:end,:); 
slp_zlb                 =   locproj_var(y,x,w,h1,H,'smooth',r,lambda);
if table == 3
    slp_zlb                 =   locproj_var_sur(y,x,w,h1,H,'smooth',r,lambda);
end
% Extract the IRF from local projection: 
irf_slp_zlb             =   slp_zlb.IR(2:H+1);
% Construct standard errors:
u                       =   slp_zlb.Y - slp_zlb.X*slp_zlb.theta;
[nwse_slp_zlb,vc_slp_zlb] = NeweyWest_slp(u,slp_zlb.X,0,0);
V22                     =   slp_zlb.B*vc_slp_zlb(1:slp_zlb.K,1:slp_zlb.K)*slp_zlb.B';
se_slp_zlb              =   sqrt(diag(V22));

%% Saving the IRFs and the corresponding standard errors. 

t                       =   0:H-1;
p                       =   0.9;
bb                      =   tinv(p,TY-P);
    
if use_data == 1
    save_LP     = [irf (irf-bb*nwse_lp(1:H)) (irf+bb*nwse_lp(1:H))];
    xlswrite('Output.xlsx',save_LP,'GDP','B2');
    save_LP_ZLB = [irf_zlb (irf_zlb-bb*nwse_lp_zlb(1:H)) (irf_zlb+bb*nwse_lp_zlb(1:H))];
    xlswrite('Output.xlsx',save_LP_ZLB,'GDP ZLB','B2');
    save_SLP     = [irf_slp (irf_slp-bb*se_slp(1:H)) (irf_slp+bb*se_slp(1:H))];
    xlswrite('Output.xlsx',save_SLP,'GDP','E2');
    save_SLP_ZLB = [irf_slp_zlb (irf_slp_zlb-bb*se_slp_zlb(1:H)) (irf_slp_zlb+bb*se_slp_zlb(1:H))];
    xlswrite('Output.xlsx',save_SLP_ZLB,'GDP ZLB','E2');
elseif use_data == 3.1
    save_LP     = [irf (irf-bb*nwse_lp(1:H)) (irf+bb*nwse_lp(1:H))];
    xlswrite('Output.xlsx',save_LP,'IP','B2');
    save_LP_ZLB = [irf_zlb (irf_zlb-bb*nwse_lp_zlb(1:H)) (irf_zlb+bb*nwse_lp_zlb(1:H))];
    xlswrite('Output.xlsx',save_LP_ZLB,'IP ZLB','B2');
    save_SLP     = [irf_slp (irf_slp-bb*se_slp(1:H)) (irf_slp+bb*se_slp(1:H))];
    xlswrite('Output.xlsx',save_SLP,'IP','E2');
    save_SLP_ZLB = [irf_slp_zlb (irf_slp_zlb-bb*se_slp_zlb(1:H)) (irf_slp_zlb+bb*se_slp_zlb(1:H))];
    xlswrite('Output.xlsx',save_SLP_ZLB,'IP ZLB','E2');
elseif use_data == 3.3
    save_LP     = [irf (irf-bb*nwse_lp(1:H)) (irf+bb*nwse_lp(1:H))];
    xlswrite('Output.xlsx',save_LP,'HOURS','B2');
    save_LP_ZLB = [irf_zlb (irf_zlb-bb*nwse_lp_zlb(1:H)) (irf_zlb+bb*nwse_lp_zlb(1:H))];
    xlswrite('Output.xlsx',save_LP_ZLB,'HOURS ZLB','B2');
    save_SLP     = [irf_slp (irf_slp-bb*se_slp(1:H)) (irf_slp+bb*se_slp(1:H))];
    xlswrite('Output.xlsx',save_SLP,'HOURS','E2');
    save_SLP_ZLB = [irf_slp_zlb (irf_slp_zlb-bb*se_slp_zlb(1:H)) (irf_slp_zlb+bb*se_slp_zlb(1:H))];
    xlswrite('Output.xlsx',save_SLP_ZLB,'HOURS ZLB','E2');
elseif use_data == 3.5
    save_LP     = [irf (irf-bb*nwse_lp(1:H)) (irf+bb*nwse_lp(1:H))];
    xlswrite('Output.xlsx',save_LP,'UNEMPLOYMENT','B2');
    save_LP_ZLB = [irf_zlb (irf_zlb-bb*nwse_lp_zlb(1:H)) (irf_zlb+bb*nwse_lp_zlb(1:H))];
    xlswrite('Output.xlsx',save_LP_ZLB,'UNEMPLOYMENT ZLB','B2');
    save_SLP     = [irf_slp (irf_slp-bb*se_slp(1:H)) (irf_slp+bb*se_slp(1:H))];
    xlswrite('Output.xlsx',save_SLP,'UNEMPLOYMENT','E2');
    save_SLP_ZLB = [irf_slp_zlb (irf_slp_zlb-bb*se_slp_zlb(1:H)) (irf_slp_zlb+bb*se_slp_zlb(1:H))];
    xlswrite('Output.xlsx',save_SLP_ZLB,'UNEMPLOYMENT ZLB','E2');    
end
 
%% Creating graph. 
 if graph == 1
        figure
        plot(t,irf,'x-b',t,irf_zlb,'*-r','Linewidth',1.5);
        h = legend('No ZLB','ZLB','Location','Best');
        set(h,'Interpreter','latex');
        hold on
        shadedplot(t,(irf+bb*nwse_lp(1:H))',(irf-bb*nwse_lp(1:H))',[0.75,0.75,0.75],'w');
        hold on
        plot(t,irf,'x-b',t,irf_zlb,'*-r',t,irf_zlb+bb*nwse_lp_zlb(1:H),'*--r',t,irf_zlb-bb*nwse_lp_zlb(1:H),'*--r','Linewidth',1.5);
        grid on
        title('Output to Productivity','Interpreter','latex','fontSize',14);
        xlabel('Horizon','Interpreter','latex','fontSize',12);
        ylabel('Percent','Interpreter','latex','fontSize',12);
        
        figure
        plot(t,irf_slp,'x-b',t,irf_slp_zlb,'*-r','Linewidth',1.5);
        h = legend('No ZLB','ZLB','Location','Best');
        set(h,'Interpreter','latex');
        hold on
        shadedplot(t,(irf_slp+bb*se_slp(1:H))',(irf_slp-bb*se_slp(1:H))',[0.75,0.75,0.75],'w');
        hold on
        plot(t,irf_slp,'x-b',t,irf_slp_zlb,'*-r',t,irf_slp_zlb+bb*se_slp_zlb(1:H),'*--r',t,irf_slp_zlb-bb*se_slp_zlb(1:H),'*--r','Linewidth',1.5);
        grid on
        title('Output to Productivity','Interpreter','latex','fontSize',14);
        xlabel('Horizon','Interpreter','latex','fontSize',12);
        ylabel('Percent','Interpreter','latex','fontSize',12);
 end

 %% Creating Table 3. 
 
if table == 1
    % Point estimates and corresponding standard errors. 
    
    Beta_h              =   [irf_slp se_slp];
    Beta_h_z            =   [irf_slp_zlb se_slp_zlb];

    Headers             =   {'Beta_h', 'Std_Errors'};
    Rowish              =   {'h_equal_0','h_equal_1','h_equal_2','h_equal_3',...
                             'h_equal_4', 'h_equal_5', 'h_equal_6','h_equal_7',...
                             'h_equal_8'};   
                         
    Table_3_Beta_h      =   array2table(Beta_h,'RowNames',Rowish,'VariableNames',Headers)
    
    Headers_z           =   {'Beta_h_z', 'Std_Errors'};
    Table_3_Beta_h      =   array2table(Beta_h_z,'RowNames',Rowish,'VariableNames',Headers_z)   

elseif table == 2
    % Pointwise test.
    
    Fslp                =   zeros(H,2);
    Flp                 =   zeros(H,2);
    chislp              =   zeros(H,2);
    
    for jx = 1:H
        Fslp(jx,1)          =   ((irf_slp(jx) - irf_slp_zlb(jx))^2/(V2(jx,jx) + V22(jx,jx)));
        Flp(jx,1)           =   ((irf(jx) - irf_zlb(jx))^2/(vc_lp(jx,jx) + vc_lp_zlb(jx,jx)));
    end
    
    chislp(:,1)         =   Fslp(:,1);
    
    for jx = 1:H
        Fslp(jx,2)          =   fcdf(Fslp(jx),1,TY-P-jx-Xp,'upper');
        Flp(jx,2)           =   fcdf(Flp(jx),1,TY-P-jx-Xp,'upper');
        chislp(jx,2)        =   chi2cdf(chislp(jx,1),1,'upper');
    end
    
    COEFF_lp            =   [irf irf_zlb];
    SE_lp               =   [nwse_lp(1:H) nwse_lp_zlb(1:H)];
    
    Headers             =   {'Test', 'P_value'};
    Rowish              =   {'h_equal_0','h_equal_1','h_equal_2','h_equal_3',...
                             'h_equal_4', 'h_equal_5', 'h_equal_6','h_equal_7',...
                             'h_equal_8'};   
    Table_3_Pointwise   =   array2table(chislp,'RowNames',Rowish,'VariableNames',Headers)
    
elseif table == 3
    % Pathwise test.
        % Construct combined variance-covariance matrix
    Vtog                =   zeros(2*H,2*H);
    Vtog(1:H,1:H)       =   V2;
    Vtog(H+1:end,H+1:end) = V22;
    R                   =   zeros(H,2*H); 
    R(:,1:H)            =   eye(H,H);
    R(:,H+1:end)        =   -eye(H,H);
    rr                  =   zeros(H,1);
    betahat             =   [irf_slp;irf_slp_zlb];
    
    for jx = 1:H
        Fslp_sur(jx,1)      =   ((irf_slp(jx) - irf_slp_zlb(jx))^2/(V2(jx,jx) + V22(jx,jx)));
        Flp_sur(jx,1)       =   ((irf(jx) - irf_zlb(jx))^2/(vc_lp(jx,jx) + vc_lp_zlb(jx,jx)));
    end
    
    for jx = 1:H
        Fslp_sur(jx,2)      =   fcdf(Fslp_sur(jx),1,TY-P-jx-Xp,'upper');
        Flp_sur(jx,2)       =   fcdf(Flp_sur(jx),1,TY-P-jx-Xp,'upper');
    end
    
    Fpath_sur            =   zeros(H,2);
    chipath_sur          =   zeros(H,2);
    for jx = 1:H
        RR                  =   R(1:jx,:);
        Fpath_sur(jx,1)     =   (RR*betahat)'*(RR*Vtog*RR')^(-1)*(RR*betahat);
    end
    chipath_sur(:,1)    = 	Fpath_sur(:,1);
    
    for jx = 1:H
        Fpath_sur(jx,2)     =   fcdf(Fpath_sur(jx,1),jx,TY-P-H-Xp,'upper');
        chipath_sur(jx,2)   =   chi2cdf(chipath_sur(jx,1),jx,'upper');
    end
            
    Headers             =   {'Test', 'P_value'};
    Rowish              =   {'h_equal_0','h_equal_1','h_equal_2','h_equal_3',...
                             'h_equal_4', 'h_equal_5', 'h_equal_6','h_equal_7',...
                             'h_equal_8'}; 
    ColumnResults 		= 	[chipath_sur(:,1) Fpath_sur(:,2)];      
    Table_3_Pathwise    =   array2table(ColumnResults,'RowNames',Rowish,'VariableNames',Headers)

end
    